{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "d3d722e4",
   "metadata": {},
   "source": [
    "## 交叉频谱和幅值平方相干性"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "be95a259",
   "metadata": {},
   "source": [
    "  创建二元时间序列。每个序列都由频率分别为 100 和 200 Hz 的两个正弦波组成。这些序列带有加性高斯白噪声，采样频率为 1 kHz。x 序列中的正弦波幅值都等于 1。y 序列中的 100 Hz 正弦波幅值为 0.5，y 序列中的 200 Hz 正弦波幅值为 0.35。y 序列中的 100 Hz 和 200 Hz 正弦波的相位滞后分别为 π/4 弧度和 π/2 弧度。您可以将 y 序列视为输入为 x 的线性系统被噪声损坏的输出。采用随机数生成器的默认设置，以获得可重现的结果。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "e75548b9",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from scipy import signal\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "fs = 1000\n",
    "t = np.arange(0,1,1/fs)\n",
    "rng = np.random.default_rng()\n",
    "\n",
    "x = np.cos(2*np.pi*100*t) + np.sin(2*np.pi*200*t) + (2*rng.random(size=t.shape)-1)\n",
    "y = 0.5*np.cos(2*np.pi*100*t - np.pi/4) + 0.35*np.sin(2*np.pi*200*t - np.pi/2) + (2*rng.random(size=t.shape)-1)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "50ee523e",
   "metadata": {},
   "source": [
    "为了防止在所有频率处获得的幅值平方相干性估计都等于 1，必须使用平均相干性估算器。Welch 交叠分段平均法 (WOSA) 和多窗谱法均适用。mscohere 可实现 WOSA 估算器。\n",
    "\n",
    "将窗长度设置为 100 个采样点。此窗长度包含 10 个周期的 100 Hz 正弦波和 20 个周期的 200 Hz 正弦波。使用默认的 Hamming 窗，重叠长度为 80 个采样点。显式输入采样率以获得单位为 Hz 的输出频率。绘制幅值平方相干性。在 100 和 200 Hz 处，幅值平方相干性大于 0.8。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "b74a701a",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEWCAYAAABrDZDcAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAABEjElEQVR4nO2dd3gc1dW437OSVtKqy7JlW5a7jRvGxsamBhOaIQkmhBBI5QuBwC+dfCTkIyGEdEIaAZKQ+gW+UAKEOGBqsOjgho17b7LcZEm2pJW0q9X9/TEz8kpW2V3tSLur8z6PHu3M3Llz75Y5c8o9R4wxKIqiKIMXz0APQFEURRlYVBAoiqIMclQQKIqiDHJUECiKogxyVBAoiqIMclQQKIqiDHJUECiuIiLrRWRBHPvbJSIXxKu/ZEBExoqIEZF0F/peICKV8e5XSS5UEKQw9k0zICIlnfa/a99Yxro9BmPMdGNMhX3dO0TkIbev2R0iMl1EXhCRGhGpE5GVInLpQI0nnojIx0VkhYg0iMh+EXlWRM4e6HEpyYEKgtRnJ3CNsyEiJwO+gRvOgPJv4EVgODAM+DJwrL8HEe8nexG5GfgV8COgFBgN3A8siud1IhhH3DUWpZ8wxuhfiv4Bu4BvA8vD9t0N3AYYYKy97wPAu1g3xb3AHZ36+TSwGzgCfMfu9wL72B3AY8DfgHpgPTC30xguABYCASAINABrwo+Htb8DeChs+1Nh176t07U9wK3Advv4Y0BxN+9FiT3nwh7er1uA/UAV8Fm7/UT7WAXwubC21wKvh23/2n7vjgErgXM6zelx4CH7+OeAAuBP9vX2AT8A0uz2afbnVA3sAL5gjyW9izEX2O/nR3uYVyaWoKiy/34FZNrHFgCVwNeBQ/Z4/qvTuXcDe4CDwO+A7E7nfhM4ADzY02cCjLXn8Rm7v2rgtrBrpQH/Y59bb7+P5faxKVhCvAbYDFw10L+vVPpTjSD1eRvIF5GpIpIGXI11QwqnEetmX4glFG4SkcsBRGQa1tPlJ4ARWDeesk7nXwY8Yp+/GLi38yCMMc9hPbE+aozJNcac0tvA7Wv/FksYjASGAKPCmnwJuBw41z5eC9zXTXdHgG3AQyJyuYiUdrrWQuC/gQuBSVjCKxqWA7OAYuDvwD9EJCvs+CIsYVAI/B/wV6AVmAjMBi7CEhAA1wMftPfPBa7s4bpnAFnAP3tocxtwuj2+U4B5WA8IDsM5/rleB9wnIkX2sZ8Ak+1zJ9ptbu90bjEwBriByD6Ts4GTgPOB20Vkqr3/Zizt9VIgH0sY+0UkB0sI/B1Lk7sauN/+fijxYKAlkf6598fxp/FvAz/Geip/EUgnTCPo4rxfAb+0X98OPBx2zIf1ZB+uEbwUdnwa0NR5DGFtH+pqjGHb7W3saz8Sdiyn07U3AueHHR+BpXGc8ORsHx+FJaS2A23Aq8Ak+9ifgZ+EtZ1MFBpBF9eqBU4Jm9OrYcdKgRbsJ2t73zXAUvv1y8CNYccuonuN4BPAgV6+B9uBS8O2LwZ22a8XAE3hfWNpBqcDgvWQMCHs2BnAzrBzA0BW2PFuPxOOawSjwo4vA662X28GFnUx/o8Br3Xa93vguwP9G0uVP7XpDQ4exLrpjcMy4XRAROZjPfnNALxY5oB/2IdHYpk8ADDG+EXkSKcuDoS99gNZIpJujGnt47g7X7ux07XHAP8UkbawfSGgVES+A3zS3vcjY8yPjDGVwBcBRKQceADr/TjDvtbKsH52RzNQEflvrKfpkVg3u3wsc5TD3rDXY4AMYL+IOPs8YW1Gdmrf01iOACW9vN8jO/Wx297X3kenc/1ALjAUS/CvDBunYJlwHA4bY5rDtrv9TMK2O39fcu3X5VhCqzNjgPkiUhe2Lx3re63EATUNDQKMMbuxnMaXAk920eTvWCadcmNMAZYd2Pnl7yfMHCMi2VgmmpiG0sW+Rjo6r4eHvd6PdXNwru3rdO29wCXGmMKwvyxjzD5jzI3GMkHlGmN+dMJAjNmLZbKY0dW1sByuEY1TRM4BvgFcBRQZYwqBoxx/D6Hj3PdiaQQlYePON8ZMj3As4bxl93V5D22qsG6m4f1V9dDeoRpLW5geNs4CY0xuWJvOn2m3n0kE19sLTOhm/yud+sw1xtwUQZ9KBKggGDxcB7zfGNPYxbE8oMYY0ywi84CPhx17HPiQiJwpIl4sM4d00UckHATGikj49241cLWIZIhIZ3v448AHReRs+9p30vE7+zvghyIyBkBEhopIl5EyIlIkIt8TkYki4rFDaj+L5UMBy6l5rYhMswXOdzt1sRq4QkR8IjIR6/10yMOy9x8G0kXkdiyNoEuMMfuBF4Cfi0i+PZ4JInJu2Fi+LCKjbFv9rT30dRTLhHaf7fvw2e/lJSJyl93sYeDb9vtTYrfvNYzXGNMG/AH4pYgMAxCRMhG5uIfTIv5MuuCPwPdFZJJYzBSRIcDTwGQR+ZQ9twwROS3Mt6D0ERUEgwRjzHZjzIpuDv8/4E4Rqce6STwWdt56LAfgI1hPqg1YNuSWGIbhmJuOiMgq+/V3sJ4Ca4HvYWkn4df+gr1vv90mfPHTr7E0mRfssb8NzO/m2gEsG/VLWJE76+w5XGtf61ks38jLWE7llzud/0u7j4PA/2I5fB2eB54DtmCZXZrpaNrpik9jmeE22PN6HMueDtbN93lgDbCKrrW4dowxP8dytH4bSxjtxTKBPWU3+QGwAngPWGv3+YNexufwTaz3420ROYb1/p3UQ/toPpPO/ALru/cC1mf0Jyw/Sj2Wn+RqLE3mAPBTLBOmEgfEGC1Mo0SOiOQCdVhO1p0DPBxXERGDNc9tAz0WRXET1QiUXhGRD9kmhxysmPK1WNE+iqKkACoIlEhYxPHFSJOwwv1UlVSUFEFNQ4qiKIMc1QgURVEGOUm3oKykpMSMHTs2pnMbGxvJycmJ74ASHJ3z4EDnPDjoy5xXrlxZbYwZ2tWxpBMEY8eOZcWK7qIge6aiooIFCxbEd0AJjs55cKBzHhz0Zc4i0u0KdTUNKYqiDHJUECiKogxyVBAoiqIMclQQKIqiDHJUECiKogxyVBAoiqIMclQQKIqiDHJUEKQwu6obWXu4r0XCFEVJdVQQpDD3V2zjvtUtaD4pRVF6QgVBCrO3ponmEFQ3BAZ6KIqiJDAqCFKYqqNNAOyp6ao6paIoioUKghSlrc2wv64ZgD01/gEejaIoiYwKghSluqGFQKgNgN1HVBAoitI9KghSlH11Te2vVSNQFKUnVBCkKFW2WSg3A/aqIFAUpQdUEKQo++qsm/9JxWlqGlIUpUdUEKQoVXXN5GWlMybfw6H6FpoCoYEekqIoCYoKghSlsraJssJshmZbH/HeWtUKFEXpGhUEKUpVnSUIhvkEgD2DxDz06PI9PL8rONDDUJSkQgVBilJ1tImRhdkM9Vkf8e5B4jB+fGUlr1aqIFCUaEi64vVK7zS2tFLnDzKyMJs8A7mZ6YMmcqjWH6RR5YCiRIVqBClIlb2GoKwoGxGhvNg3aNYS1PkDNASNJtpTlChQQZCCVDqCoDALgNHF2ew+kvr5howx1PmDtLZBc7BtoIejKEmDCoIUxNEIRhZmAzBmSA57a5toa0vtp+T6llZa7TnW+jXjqqJEigqCFKSqrol0jzAsz9IIyot9BFrbOFjfPMAjc5faxuM3/zq/OgoUJVJUEKQg+2qbGF6QRZrHCh0dU+wDUj+EtDbs5l/XpBqBokSKCoIUpKqumTLbLAQw2hEEKe4wDjcHqUagKJGjgiAF2WcvJnMoK8rGI6kvCOpUEChKTKggSDFaQ20cONbc7igGyEjzMLIwO+UFQU2jmoYUJRZ0QVmKcai+hVCboawou8P+0cW+lM9CWucPIAJpAkdVI1CUiFGNIMXY1yl01GHMEF/Kry6u9QcozM4gL0M0fFRRokAFQYpR1WkxmUN5sY8jjQEaWloHYlj9Qq0/SJHPS06G+ggUJRpUEKQYlbXdaATFOUBqh5DWNgYo9GWQkyHUNakgUJRIUUGQYlTVNVHky8Dn7ej+GQwhpLX+IMU5XnK90iGCSFGUnnFVEIjIQhHZLCLbROTWLo6PFpGlIvKuiLwnIpe6OZ7BQFVd0wmOYoDRQxxBkLo5h+r8AQp9XksjUNOQokSMa4JARNKA+4BLgGnANSIyrVOzbwOPGWNmA1cD97s1nsHCvromRhacKAgKsjMoyM5IcY0gQJEvg1zbNKQZSBUlMtzUCOYB24wxO4wxAeARYFGnNgbIt18XAFUujiflMcawr7bpBP+AQyqHkDYFQjQH22yNAAKtbZqBVFEixE1BUAbsDduutPeFcwfwSRGpBJYAX3JxPCnPsaZWGgMhRnVhGgLLPJSqIaROuGhxjpfcDOmwT1GUnhnoBWXXAH81xvxcRM4AHhSRGcaYDo9yInIDcANAaWkpFRUVMV2soaEh5nOTgT3HQgDU7ttBRcUeoNOc6wPsrQny8tKleEQGaJTusNue+74dW0gLtQDCS6++yej8tIEdWD+R6t/trtA5xw83BcE+oDxse5S9L5zrgIUAxpi3RCQLKAEOhTcyxjwAPAAwd+5cs2DBgpgGVFFRQaznJgMvbTgIb67gorPmckp5IdBxzgd8e3hm51omnTKfcjuKKFV4Y1s1vPkOZ8+bzbvvrgaamTj9FM6cUDLQQ+sXUv273RU65/jhpmloOTBJRMaJiBfLGby4U5s9wPkAIjIVyAIOuzimlKa7VcUOTghpKpqHHDNQkc8KHwVdVKYokeKaIDDGtAJfBJ4HNmJFB60XkTtF5DK72deB60VkDfAwcK3RUI+YqaprwpvuYUiOt8vjTgjp7lQUBI2OIMggN8Pap4JAUSLDVR+BMWYJlhM4fN/tYa83AGe5OYbBhJN+2uPp2v4/oiCbdI+kZAipU5TGWUcAmoFUUSJFVxanEPvqmhjZKcdQOGkeYVRRdkqmmaj1B8jNTMeb7sGbJmSmezQDqaJEiAqCFKKqm8Vk4YwekpOSGkGdP0ihL6N9u8jn1fBRRYkQFQQpQqC1jUP1LV2mlwhndHFqFqipaQxQ5DvuGyn0ZaiPQFEiRAVBinDgaDPGdB8x5DCmOIejTcGUM5vU+QMUhTnJC7IzNAOpokSICoIUwQkdHdWLIChP0SykVi2C46YhSyNQ05CiRIIKghShtzUEDs5agt0ploXUSjh3XCMo8nnVNKQoEaKCIEVwKpMNL+g+agjC01GnjkYQDLVR39zawVlc4MvQDKSKEiEqCFKEfbVNDM3LJCuj59w6uZnpDMnxptTqYufJvzjMR1CY7dUMpIoSISoIUoSqo92nn+7M6CGplY7a8QUUdjANWdqBhpAqSu+oIEgRrFXFPZuFHEYX+1LKNOSsKu7sLAZNM6EokaCCIAUwxlglKiPUCMqLfFTVNdHWlhr28/CEcw4F2dZrTTOhKL2jgiAFqGkM0Bxsi9g0lJ+dTpuBpmDI5ZH1D+0J53I6LigD1QgUJRJUEKQATuhopBqBz2vlGmwMtLo2pv6kK9OQox2oIFCU3lFBkAJURbiGwCEn04os8rekhkZQ5w/gTfeQHRYx1a4RqGlIUXpFBUEKsK+uGRjMGkGAIl8GElZ+MysjTTOQKkqEqCBIAfbVNuHzpnVYUNUTPq+tEQRSQyOoaQx2cBQ7FPoyNHxUUSJABUEKcPBYM8MLsjo8EfdEu0bQkhoaQV2n9BIOmmZCUSJDBUEK0NDSSl5m5MXm2n0EKaIR1PoDFOWcqA1pBlJFiQwVBClAUyBEtrfn1BLh5KScRhDssKrYodCXoT4CRYkAFQQpQFMw1CFipjdSyUfQ1maoa+qYgtpBq5QpSmSoIEgB/IHWdrt/JORkpk7UUH1zK6E206WPQDOQKkpkqCBIAaI1DWWme/BIaqwj6Cq9hINmIFWUyFBBkAL4g6F2c08kiAg53vSU0AjaBUEXzuJCzUCqKBGhgiAFaApE5yMA8GWmpYRG4ISHduUsLtJ8Q4oSESoIkpxQm6GltS0q0xCQMhpBjZ1wrrgrH4FmIFWUiFBBkOQ4GUSjMQ2BrRGkQNRQjz4CWyPQEFJF6RkVBElOk30zj9o05E1PiXUEdf4gHoG8rBOjpo77CFQQKEpPqCBIctoFQRThowA53tTRCAp9XjyeE9NrtKeiVtOQovSICoIkxx+0nuqjNw2lho/AyTzaFZqBVFEiQwVBkhOraSjHmxpRQ7XdZB510AykitI7KgiSnOOmoRh8BCmiEXQVOuqgGUgVpXdUECQ5jp0/WtNQjh01lOzpF+r8XecZctAMpIrSOyoIkhx/rOGj3vT2NQjJijGGWn+A4pyeTUPqI1CUnlFBkOQ02xpBVgw+AkjuDKRNwRAtrW09moYKszUDqaL0hquCQEQWishmEdkmIrd20+YqEdkgIutF5O9ujicV8QecqKHowkd9mclfk8BZH9CTaagwRzOQKkpvRHf3iAIRSQPuAy4EKoHlIrLYGLMhrM0k4FvAWcaYWhEZ5tZ4UpVYTUNOcZpk1ghq7fQSvWkETgbSaB3qijJYcFMjmAdsM8bsMMYEgEeARZ3aXA/cZ4ypBTDGHHJxPClJcyCEiJVaOhp8drnKZI4ccqKBevMRgGYgVZSecE0jAMqAvWHblcD8Tm0mA4jIG0AacIcx5rnOHYnIDcANAKWlpVRUVMQ0oIaGhpjPTVQ272jB64FXXnmly+PdzXlLjaUJvLVsFcd2JOeT8tv7LSG2df1q/LuPC8LwOVcesNq89OqbjM5PznlGQip+t3tD5xw/3BQEkV5/ErAAGAW8KiInG2PqwhsZYx4AHgCYO3euWbBgQUwXq6ioINZzE5UXateSV32g23l1N+eSfUdh2etMnDqdBdOHuztIl9jz1i5Ys56Lzj2LoXmZ7fvD55y5/Qj3rX6bidNP4cwJJQMz0H4gFb/bvaFzjh9umob2AeVh26PsfeFUAouNMUFjzE5gC5ZgUCKkKRCKOmIIjper9Cexaai20alF0IOzWDOQKkqvuCkIlgOTRGSciHiBq4HFndo8haUNICIlWKaiHS6OKeVoCkRXnczBCR9tTOI0E7X+AHmZ6WSkdf811gykitI7EQkCESkVkT+JyLP29jQRua6nc4wxrcAXgeeBjcBjxpj1InKniFxmN3seOCIiG4ClwC3GmCOxTmYw4g+Gos48CsfDR5NaI/AHKOrBUQyagVRRIiHSO8hfgb8At9nbW4BHgT/1dJIxZgmwpNO+28NeG+Bm+0+JgaZAK74YTENOkrrk1gh6Ti8BmoFUUSIhUtNQiTHmMaAN2p/2k/cOkkI0BUMxxceneYTsjLSk1gjqekk456AZSBWlZyIVBI0iMgQwACJyOnDUtVEpEeMPxCYIwEo815jMC8p6qEUQTmG2ZiBVlJ6I1DR0M5ajd4Id8z8UuNK1USkR0xQIxWQaAisthT+JU0zUNQZ79RGApRFoBlJF6Z6IBIExZpWInAucBAiw2Rijv6wEIFbTEFhpKZJVIwi0tlHf0tpjURqHQl8Gu6r9/TAqRUlOIo0a+gKQa4xZb4xZB+SKyP9zd2hKJPTNNJSetD4CJwooUtOQ+ggUpXsi9RFcH77a184NdL0rI1IiJtRmCLS24cuIbYG4z5uWtFFDjs0/UmexZiBVlO6JVBCkiYg4G3Zm0d5/gYqrOE/z2d7Y1gXmeJNXI3Ayj/aUcM6h0Hc8A6miKCcS6R3kOeBRETlfRM4HHrb3KQNIU9CpVxyjRpCZvBqBY+rpKb2Eg2YgVZSeifQO8k3g88BN9vaLwB9dGZESMU7h+lijhpJaI2gvShOBRpBtCYI6f5CRhdmujktRkpFIo4bagN/af0qC4BSViTlqKInXEThP95FFDWmaCUXpiYgEgYicBdwBjLHPEawMEePdG5rSG8dNQ7FrBIHWNoKhth4TtyUidf4gWRmeiOauGUgVpWciNQ39CfgasBJNLZEw9NU05AsrYF+QnVyCoKYxEJE2AJqBVFF6I1JBcNQY86yrI1GixjENRVu43iG8JkFBdu9O10Qi0jxDYK0jADUNKUp3RHoHWSoiPwOeBFqcncaYVa6MSomI46ah2J7mfUlckyCSzKMO2V7NQKooPRGpIHBqDc8N22eA98d3OEo0NLWvI4hRI/Amb02Co01BJg3Ljbi9ZiBVlO6JNGroPLcHokSPv68+gszk1Qjqm4PkZUUuADUDqaJ0j2sVyhT36Wv4aDJrBPXNreRlRe7XKNAMpIrSLZEal/+KVVZypL29BfiqC+NRoqA5GEIEMtNjTDHhaARJtpagNdSGPxCKSiMo8mWoj0BRukErlCUxfrsWQVgaqKhwoo2SrSaBY8qKRiMozvFypLGl94aKMgjRCmVJTF9SUMNx01CyaQTHmq0n+2g0gmF5WRxpDBAMaeI5RemMVihLYpr7UJQGjvsWkk0jqG+2xpsfhSAozc/CGKhuaGFEgeYbUpRwev0l2Smnz7X/tEJZAuEPtMZciwDAm+7Bm+ZJOo2gvl0jiNw0VJqfCcDBYyoIFKUzvZqGjDEh4BpjTKtToUyFQGLQV9MQWCGkyRY15GgE0ZqGAA4ea3ZlTIqSzET6S3pDRO4FHgUanZ26snhgaQ6GyI5xDYFDjjc96dYR1LfErhEcUkGgKCcQqSCYZf+/M2yfriweYPyBEMPz+5YjyOcdHBrBkNxMPGKZhhRF6YiuLE5imuJiGkpPQh+BJQhyMyMXBGkeYWheppqGFKULdGVxEuMPxMM0lJZ0UUPHmoN40zxkRTn30vwsDtarRqAondGVxUlMUzDUnkE0Vnze5NQIojELOQzLy1IfgaJ0ga4sTmIs01Ds4aNgpZlIRh9BLIKgNF9NQ4rSFbqyOElpDbURCLX12TTkS8aooeZgVBFDDqX5WdT6g7S0Jtd8FcVtdGVxkuIUpemraSjHm9Ze1yBZ6ItGAHC4voVRRb54D0tRkpZIo4ZWiYiuLE4gmvqYgtrBl5mOPxiirc3g8cSWvK6/qW8OUlKSE/V5w/KdRWUqCBQlnGgeq+YBY+1zThURjDF/c2VUSq8cr1fcd43AGGhuDcVc+7i/aYiyFoFDqb26WB3GitKRSMNHHwTuBs4GTrP/5vZ4knXeQhHZLCLbROTWHtp9RESMiPTap2LRXpSmrz4COxY/mfwEfTUNqcNYUToS6a9pLjDNGGMi7dhOVncfcCFQCSwXkcXGmA2d2uUBXwHeibRvJbxwfd81AnCqlGX2dViu09ZmaAjEphEU+bxkpImuJVCUTkQaNbQOGB5l3/OAbcaYHcaYAPAIsKiLdt8HfgroY1oUNLWbhvpmznHOTxaNoCHQijHRpaB28HiEYXlZqhEoSid6/DWJyL+xQkbzgA0isgxof5wyxlzWw+llwN6w7Upgfqf+TwXKjTHPiMgtPYzjBuAGgNLSUioqKnoadrc0NDTEfG6iseqgFemz4b138e/uXivobc7bqi0B8MY7yzlU1Dftoj840mQVltm3ezsVoT1dtulpzlmmhc2791NRUefSCAeGVPpuR4rOOX709lh1d9yvaCMiHuAXwLW9tTXGPAA8ADB37lyzYMGCmK5ZUVFBrOcmGkdX74N3V3P2GfOZOCy323a9zTlvdy2seJPJ02dy7uShLow0vmw6cAxeeY3TTpnBgpNHdNmmpzk/sncl2w83sGDBuS6Osv9Jpe92pOic40ePgsAY84rzWkRKsZzEAMuMMYd66XsfUB62Pcre55AHzAAq7Jq7w4HFInKZMWZFZMMfvDTFK2ooM7mqlMWScC6c0vxM3txeHc8hKUrSE2nU0FXAMuCjwFXAOyLS24Ky5cAkERknIl7gaqxFaQAYY44aY0qMMWONMWOBtwEVAhESr6ihZKtbXB9DveJwhuVncay5tV2QKooSedTQbcBpjhYgIkOBl4DHuzvBGNMqIl/ESlaXBvzZGLNeRO4EVhhjFnd3rtI78Yoa8nWIGkp8jtciiK0OQ6m9qOxQfTNjhkS/KE1RUpFIBYGnkynoCJGVuVwCLOm07/Zu2i6IcCwKlmnII5CZHmngV9fkJNk6gmMxFK4PJ7x2sQoCRbGI9Nf0nIg8Dzxsb3+MTjd4pX9xahHY/pWYyUz34JFk0giiL1MZTmm+1i5WlM70Fj46ESg1xtwiIldgrSwGeAv4P7cHp3RPU7C1zymoAUQkqeoW1ze3ku4RsjJi04RKtYi9opxAb3eSXwHfAjDGPAk8CSAiJ9vHPuTi2JQeaAr0vSiNgy+JahJYKajTY9aE8rPTyUz3cEhXFytKO709VpUaY9Z23mnvG+vKiJSI8MdREOQkUZWyWBPOOYiIVbJSNQJFaac3QVDYw7HsOI5DiZKmYCjqmr3d4ctMnrrFsSacC2eYFrFXlA70JghWiMj1nXeKyOeAle4MSYmEuJqGvOk0Jo1pqO+CoDQ/i0PH1DSkKA69/aK+CvxTRD7B8Rv/XMALfNjFcSm94A+EKPTFbiIJJ8ebxpHGQFz6cptjzUHKi/tWVGZYfiYVm1UjUBSH3lJMHATOFJHzsNJBADxjjHnZ9ZEpPRJf01A6e2r8cenLbeKlETQGQjS0tMacqkJRUolIS1UuBZa6PBYlCuJqGspIa09ZkejUNwfJ6+PN21lUduhYM7lDu0/YpyiDhb4tS1UGDH+gNW6lJXMy02lMAmexMYaGlr5FDUH4WgL1EygKqCBIWuJqGvJaGkEUBegGhMZAiDYTe8I5h2Fh+YYURVFBkJQEQ20EQyZ+6wgy02ltMwRCbXHpzy36ml7CQWsXK0pHVBAkIU7m0fiFjzo1CRLbT3A882jfNILczHR83jQ1DSmKjQqCJMTJpd/XFNQOx2sSJLafoK+1CBx0dbGidEQFQRISr6I0Dj6nSlmCRw4d62MtgnCG5WXqojJFsVFBkITEq0ylQ7tGkOCRQ/V9rEUQTml+FgfVWawogAqCpKQpaN0Q45GGGsKrlCW2RhAvZzFYDuODx5oTPlJKUfoDFQRJSLxNQ8erlCW2RtAQJ2cxWBpBc7Ct3dykKIMZFQRJSLxNQ8mjEbSS5pG4zLt9LYE6jBVFBUEyEq/C9Q7tGkESRA3lZsZelCac0rzjtYsVZbCjgiAJiXvUUBKtI4iHWQi0drGihKOCIAnxx900lBwawbHm+GULHeasLtbIIUVRQZCMNMfZNJRmF4NPfB9BkPw4RAyBJfzystJ1LYELNLS08vtXthNq04isZEEFQRLiD7TiEfCmxe/jy/EmfgbSeJqGAF1d7BJL3tvPj5/dxKo9tQM9FCVCVBAkIVbh+vg4TR18mYlfk6C+JRhnQaC1i91g88F6ACprk6PYkaKCIClpDobiZhZySB6NID6mIbDqEmjUUPzZYguCvTVNAzwSJVJUECQh/jhWJ3NwahIkKsaYuJuGhuVncaheVxfHmy2qESQdKgiSEH8gFLfQUYeczPSEjhpqCoYItZn4agT5mQRDhlp/MG59DnaONgXbtSzVCJIHFQRJSFMg/qYhnzctodcRxKsWQTjD8nQtQbzZamsDhb4M9qpGkDSoIEhCmoLxNw3leBNbI4hXLYJw2ovY16ufIF44juLzThrG/qPNtCZ41TvFQgVBEuKGaSjRo4aOp6COp2lINYJ4s/VgAzneNOaNKybUZth/VN/bZEAFQRLSFGiNWwpqh0SPGnLDNDTUzjekiefix+YD9UwqzaO8yAdAZa36CZIBFQRJSFMwhC/eGoE3nZbWtoRV5evjWJ3MISsjjUJfhoaQxpGth+qZXJpLeXE2gPoJkgRXBYGILBSRzSKyTURu7eL4zSKyQUTeE5H/iMgYN8eTKvhdcBbnOOUqg4lpHnJ8BLlx1AjAWUugGkE8ONLQQnVDgMmleYwoyEYEKmtUECQDrgkCEUkD7gMuAaYB14jItE7N3gXmGmNmAo8Dd7k1nlTCnagh6wabqJFDbpiGwEo+d1CdxXFhy8EGACaX5uFN9zAiP0tNQ0mCmxrBPGCbMWaHMSYAPAIsCm9gjFlqjHEeGd4GRrk4npQgGGqjtc3E3TTkaASJGjlU3xxEBHLj7Bspzc9SH0Gc2HrIihiaXJoHwKhin5qGkoT4/qo6UgbsDduuBOb30P464NmuDojIDcANAKWlpVRUVMQ0oIaGhpjPTRQag9Yq2Kq9u6io2Ndr+0jnvOOQJQBee/Md9hbEV8jEg43bW8hKg1dffaXXttF8zsG6AAeOBnn+paVkpscvd1N/kwjf7Yr1LfjSYeOqt9gkQkZLCxtrQq6NKxHm3N+4NWc3BUHEiMgngbnAuV0dN8Y8ADwAMHfuXLNgwYKYrlNRUUGs5yYKB442w3/+w4ypk1kwv3eXSqRz9m6rhlXvMOXkWZw+fkgcRhpf/n1oDUV11RHNJZrPuXXYQf61fQXFE0/htLHFfRvkAJII3+37N73FtDLDeeedCcC7wS28+fJWzjj7HDLT4/9wkQhz7m/cmrObpqF9QHnY9ih7XwdE5ALgNuAyY4waa3vBb5tu4p5ryC744k9g01A8I4YcZpYXALBmb13c+x5MGGPYcsgKHXUYVZSNMVBVp6a3RMdNQbAcmCQi40TEC1wNLA5vICKzgd9jCYFDLo4lZWivV5wR73UEto8ggZ3F8XYUg5VmYmRBFmsqj8a978HE4foW6vxBTirNbd9XXuysJVA/QaLjmiAwxrQCXwSeBzYCjxlj1ovInSJymd3sZ0Au8A8RWS0ii7vpTrFpCsS3OplDwmsEca5FEM4p5YWqEfSR8IghB0cQaPK5xMdVH4ExZgmwpNO+28NeX+Dm9VOReNcrdkgGjWB8SW7vDWNg5qhCnl13gNrGAEU5Xleukeo4qacnDz8uCIbnZ5HuEdUIkgBdWZxkHDcNubSOIFE1ApdMQwCnOH6CyjpX+h8MbDlYT3GOl5LczPZ9aR5hZGE2e3UtQcKjgiDJcMs05E33kJEmNCZg4jmrKI07zmKAk8sKEIH31E8QM1sO1jNp2Ika26iibPbq6uKERwVBkuGWacjqMx1/Aiaea2ltIxgyrmkEeVkZTBiaq36CGDHGsPVgAyeFmYUcyot8uro4CVBBkGQ4piFfnKOGwPITJKJGcDwFtXsurVNGFbKmsk7LVsbA/qPN1Le0dggddSgvzqa6oaVdk1USExUESUaTbcOPt2kIrMihRPQRuJVwLpxTyguobghQpfnzo6bdUdylaciKHNpXp+ahREYFQZLhD4RI8wgZafFPh5DjTUvIqKH2hHOZ7vgIwNIIQBeWxUK7IOhGIwANIU10VBAkGf6AVYtAJP6CwOdNVI3Ancyj4UwZkUdGmmjkUAxsOdjA0LzMLkNvHY1Ak88lNioIkozmYPxTUDvkZCaqRuDUK3ZPI8hMT2PaiHzVCGJgy0GrGE1XDM3NxJvuUYdxgqOCIMlwoyiNw2DWCMBaYbxu3zFCbeowjpS2NitiqCuzEIDHIxpCmgSoIEgy3Chc75CTmZhRQ8dsjSCeheu7YuaoQhpaWtlxuMHV66QS++qaaAqGuhUEYJmHVCNIbFQQJBnNwZArawggcdcROBqBm1FDALPsFcar1TwUMZsPdO8odigvylYfQYKjgiDJ8Ada3fMReNPwB0O0JZhppL65lRxvGmkedwvHjC/JJTczvd9WGDe0tFLbGKCxpZXWUFu/XDPebLGrkk3qxkcAVvK5On+w3dejJB4JUZhGiRx/IERxTmbvDWPAl5mOMdDcGmrPPZQIuJleIhyPRzi5rKBfIof2HPFz0a9eoTl4XAB4xHJaZ3vTuGpuOV+/aDIZaYn9rLb1YAMjCrJ6NNuNKrJCSCtrm5g6wv3PUYmexPm1KxHhpmnI+TFX1wcYPSRxvhpuJpzrzMzyAv78+k5aWkOuVNVyuHfpVoyB73xwGqG2NlqCbbS0thEItbG3xs/vXtnOOzuP8JtrZreHYMYDYwy/fGkrR/0B7rhsep/DkDcfqO/RLARWmgmAvTV+po7I79P1FHdInF+7EhFuOovnjCkC4O0dRxg9JH43n77iZi2CzswaVUgwZNi4v55Z5YWuXGPPET9PrNrHp88Yw3Vnj+uyzdPvVfGtJ9Zy6a9f42cfPYWLpw/v83Xb2gy3PbWOh5ftAawoqStOHRVzf6E2w7bDDZw1sefSpuEagZKYJLbeqZxAk4vho5NLcxmWl8mrWw+70n+sNDS39otpCKybI8B7LpqH7l26lXSPcNO5E7pt88GZI3n6y2czZkgOn39wJd/793paWmOP6GoNtfHf/1jDw8v2cNOCCcwZU8SdT2/gcH3s1WF3H2kk0NrWZY6hcIpzvPi8aeowTmBUECQRxhj8LpqGRISzJ5XwxrbqhHIY96dpaERBFiW5ma5FDjnawMfnj2ZYflaPbccMyeHxm87gs2eN4y9v7OIjv30zpiIvwVAbX3lkNU++u4+vXziZby6cwk8/cjL+lhDfXbwu1qm0VyU7qRdBICKUF/k0zUQCo4IgiQiGDKE245ogAHjfpKHU+oOsrzrm2jWi5Vg/CgIRYVZ5gWsrjO9buq1XbSCczPQ0bv/QNP7w6bnsqvbzg6c3RnW95mCImx5ayTNr93PbpVP50vmTAJg4LI+vXDCJJWsP8Ny6/VHPA+DVrYfJzkjr1UcAlnlIK5UlLioIkggnlW+WSz4CgLMmlgAklHmov6KGHGaOKmRHdWP7QrZ4YWkDlVwzr3dtoDMXTivlo3NH8fKmQxxtimxcTYEQ1/9tBS9tPMT3F03n+veN73D8hveNZ9qIfL791Hrq/IGoxhMMtfHs2v1cMK00IlNlebG1qEzTfMeOm+G3KgiSCH/QWljlZmjn0LxMpo7I5/Wt1a5dIxoCrVY0TV5m/8U1nFJeiDGwLs7rCe5bug2PR7hpQWTaQGcun1VGwL4BR8J/P76G17dVc9dHZvKpM8aecDwjzcNdV86k1h/g+1FqGq9vrabWH+SyU0ZG1H5UUTYNLa0RCzGlI4eONfO+u5bySqU7758KgiSiycXqZOGcM6mEFbtrEiLv0PGEc/0nCGaWOTWM4ycIHG3g4/NGUxqlNtA+rlEFjC/J4anV+3ptu7fGz5K1+7nx3AlcdVp5t+1mlBVw47njeWJVJRWbD0U8lsVrqijIzuDcyUMjat+ehVT9BDFx+7/W0xgIcVKRO799FQRJhL8fTENgCYJgyPDOzpq49muM4bl1B6KyFR9PONd/pqGiHC9jhvii8hNU1vq58cGVPLmqkmAXq4T7qg2A5b9YNKuMt3fUUFXX8w31oXd24xHh02eM6bXfL71/EhOG5nDbP9fREEGKkaZAiBfWH+CSGcPxpkd2C2mvS6B+gqh5bt1+nlt/gK9eMInhOe7cslUQJBHtZSpd1ghOG1tMZrqH17bE1zz0woaD3PjQSi78xav86fWdEWX57K/Mo52ZOaowqhDSHz+7iefWH+Dmx9aw4GcV/O+bu9o1uHhoAw6Xz7ZMMYvXVHXbpjkY4rHle7lg6jBGFGT32mdWRhp3XXkKVUeb+Omzm3pt/59NB2kMhCI2C8FxjUAdxtFx1B/kO/9az7QR+Vx/zvjeT4gRFQRJhJuF68PJykhj3rhiXt8WP4dxUyDEnf/ewOTSXE4fX8z3n97AFb99k00Heo5O6o9aBF0xb2wRVUebeS0Cp/nqvXU8895+vvz+ifz52rmMKMjiu4vXc9ZPX+ae/2zl7hc24/EIN0YYKdQTY4bkcOroQp56t3vz0NPv7afWH+TTXfgFumPOmCI+c8ZYHnpnd6+fyeLVVQzLy2T++J4XkoVTkJ1Bfla6moai5MfPbqSmMcBdV850Nd2ICoIkoj+ihhzOmVTCloMNHOihhu/h+hb+559rI8o1f9/Sbeyra+IHl5/Mn689jV9fPYu9NX4+eM/r3P38ZpqDXS+WOjZAGsFH55YzriSH7zy1rtuxgWXu+vGSjZTkernh3Am8f0opj990Jv+48QxmlRfyixe3sHhNFR+fN5rhBX3TBhwun13GpgP1bNzf9Q37wbd2MWFoDmdOiPxGDfCV8yeRm5nOz57b3G2bo01BKjYf5oMzR0adBNBKR538GkFTIMTayqM8sbKSHz+7kc/+dTln//RlZt7xPF9++F2eXbu//bfaF97cXs0jy/fyuXPGMcP2W7mFpphIIprao4bcFwRnTxwKbOL1bdVcOafrNAQ/e34Tj62oZMWuGp646cxun9p3HG7ggVd3cMXsMuaNKwZg0awyzpk0lB88s4F7l25jydr93PvxU5k2smMumvp+qkXQmayMNL6/aAaf/NM73F+xnZsvnNxlu6WbD/HOzhruXDSd3LDIptPGFnPatcVsOnCMJWsP8F9njo3b2D5w8gju/PcGnlq974TcPWv21rGm8ih3fGha1HmEinK83LRgAnc9t5llO2vaP6twnl9/gECojctmRW4Wcigvzmb74caoz0sk/vz6Tn7wzAYcq2ZGmjBhaC6zRxfhTfPw8qaDLF5TRVaGh/NOGsbCGcN5/5RhUWu0TYEQ33pyLWOG+PjaBV1/9+KJCoIk4rhpyP2PbcrwPEpyM3lt6+EuBcGGqmP8Y2Ul50wq4c3tR/jao6v5/afmnvCUaIzhu4vXk5nu4dZLp3Q4Vpzj5RdXzeLyWWXc8vgarv/bCpZ8+RwKfMd/NAPlIwA4e1IJl88aye8qtrNo1kgmDO2YajnUZvjJs5sYO8THNfNGd9nHlOH5TBke30RrQ3Ized/koSxeXcU3L56CJ+w9f/Dt3fi8aVzRjfDujf86cxz/++YufvLsRp646cwThMni1VWMGeLjlFHRP6GWF/l4ZcthjDHt/baG2lhTWceq3XV8ZM4oiruoe5wo7Dni56fPbeKMCUP4xPwxTC7NZcyQnA4mm9ZQG8t21vDsugM8t/4Az647QHZGGt9YeBKfOWNsh8+qJ3710hZ2H/Hz9+vn94sFQE1DSYSjbrqVaygcj0c4e+KQLtNNGGP40ZKNFGRncO81p3L7B6fx0sZD3P3CiSaF59Yd4LWt1dx80WSG5XVtGnnf5KH87pNzOHismW88sabDoqP+KkrTHbd9YBpZGR5u++faExZDPbGyki0HG/jGwin9ni768tll7D/a3CGyq7YxwL/XVHH57LKYNahsbxpfu2Ayq/bU8cKGgx2OHapv5s3t1Xxo5siYspaOKsqmOdjGhv3H+MeKvXzh76uY84OX+Mhv3+KHSzZy44Mru4y4SgSMMdzx7/Wke4Sff3QWl548gonD8k743NPTPJw5sYTvXz6Dd751Po/feAbzxxfzvX9v4Jo/vM2eI72bxtZWHuUPr+3g6tPKOXNCiVtT6oAKgm5YuukQ77trKefdXcGH73+D//rLMm5+dDV3/nsD91ds49Cx7m3nbtEuCPrhCQHgnElDqW4IsLGT87Bi82Fe31bNV86fRIEvg0+fMYaPzx/Nbyu28893K9vb+QOt3Pn0BqaOyOdTp/ccxjh7dBG3XjKF59cf5G9v7W7fX98cJDsjbcDy8g/Ny+Sbl0zh7R01PLnquIO2KRDiFy9uYVZ5IZfM6Htm0Gi5cGopOd40/hW2puCxFXtpaW2LKGS0J66cM4oJQ3O467lNHQrmLHlvP20GFsVgFgJrdTHAB+55nVsef4/lO2u4aFop9358Nj+54mSW7arhh89Et7Ctv3hxw0Fe3nSIr14wOWJfj8cjzB1bzF+uPY27PjKT9VXHWPjrV3no7d1drrAOhtpYubuGWx5fQ0luJt+6dGq8p9EtahrqgqWbDvH5B1cytsTH5NI8jjYFOdzQwpaDDRxtCtLQ0sr9S7fz1Qsm8Zkzx/bbTcofDJHukYhjt/vK2ZOsp5HXtlYzfaRlCmgNtfHDJRsZV5LDJ+ZbNxwR4XuXTWf7oQa++cRaxg7JYfboIn7z8jb2H23mN9fMJj2C9+i6s8fx1vYj/PCZjcwZU8SMsgIaWloHTBtwuOa00TyxspIfLtnI+6cMoyjHy5/f2MmBY838+upZfc7pHwvZ3jQunjGcZ9bu547LptNmDA8t2828scV9NkWlp3m45eIp3PjQSh5fWcnVttlr8ZoqpgzP6zXbaHfMG1fMFbPLmDAslwUnDWXaiPwO7932ww384bWdzCgr6NYvFSuhNoNAxKaZcJoCIb5nR7xde9bYqM8XEa46rZyzJpXwzcff49tPreO5dQf48RUnU+sP8Ob2I7y1/QjLd9XgD1i/8d99cg4F2f3nF1NB0AlHCJw0PI+HrpvfwV7tsLO6ke/9ez0/eGYjj63Yy52LZnB6FKF0seJmCuquKM3P4qTSPF7fWt0e+vjw8r1sO9TA7z81p4NAykjz8NtPzmHRfa9zw4Mr+dXHZvHH13Zw5ZxRzB17otOxK0SEuz96Cpfe8xpf+Psqnv7S2f2aebQ7PB7hhx8+mQ/+5nV+8uwmvnnJFH5XsZ3zpwyLKoQy3lw+q4wnV+2jYvMhNleH2FvTwjcuntL7iRFw8fRSZo8u5JcvbWHRrDKqG1pYtaeObyw8KeY+87Iy+MXHZnV7/JsLp7C+6hj/88+1nFSax8kx+CG6Ymd1I5/73+X4vOk88Ok5Ea2tCOfepVvZV9fEozec3qeHvrLCbB68bh5/X7aHHz6zkXPuWtp+bOKwXK6cM4ozxg9h/vgh/e4rUdNQGJEIAYBxJTn85drTeOBTc/AHQlz9wNt8+eF3OeiyuajJxaI03XH2pBKW7aqhORiivjnIr17cwrxxxVw0rfSEtsU5Xv70mdPwt7TyyT+9Q1ZGGrdeEt2NqSjHyz3XzKaytolvPbmWY/2ccK47po7I53Nnj+PRFXv50sOraAy08s0o5xZvzpwwhKF5mfzz3X38Z08rJbmZcSlgA5ZQvnXhFA4ea+Gvb+5qX8D2oZmxmYUiIT3Nw2+umc3Q3Ew+/+AKjjTEXivBYdnOGj58/xvUNAbYcbiBy+97g7VRpA7Zduh4xFs8hL6I8In5Y3j+q+/jS++fyK+vnsWy/zmfl24+lzsXzeCSk0cMiMN80AiCXdWNrDnc2m0Gv6WbLSEweXhuj0LAQUS4aPpwXrr5XL58/iSeW3+A999dwRf+voq/vrGT9VVHI1o525mj/iD/Wr2PF9YfOCEW2c1aBN1xzqQSAq1WJMT9Fds50hjg2x+Y2q05ZHJpHvdcM5t0j/CtS6ZSkht9feXTxhZz84WTefq9/byzs4b8AdYIHL5ywSTKCrN5Y9sRPjqnPKL0y26SnubhQzNH8vKmQ6w9HOLj88rjajacP34I508Zxv0V23h8ZSVzxhS12/ndYkhuJr//1ByONAb4wt9XdfBRGGPYuP8YP39hMx/8zWv84b0WdlZ3H4765KpKPvHHtynO8fLUF87iif93JukeDx/9/ZsRpd62It7WkZWRFnd7fXmxj69fdBKLZpVFnYnWDVz9hYnIQuDXQBrwR2PMTzodzwT+BswBjgAfM8bscmMsT63ex69WtvDrVS8wfWQB88cVM3/8EOaNLebdvbV8/sGVTCqNTAiEk5WRxs0XTuYjp5Zxz3+28eb2ap55z/qS5WWmc+qYIk4bW8TUEfmMK8mhvNh3gnp5pKGFFzYc5Nl1B3hzWzWttgDJyvBwzqShXDStlPOnltIUaCW7n4vKzx83BG+ah0eX7+XFjQf58OwyZo4q7PGc86eWsua7F/UpzPWmcyfw9o4jvLa1esBNQw4+bzp3XTmTn7+wma91s66gv/nw7DL+/MZOPAIfn983J3FXfGPhFBb++lXqm1u5No5rIXpiRlkBP77iZG5+bA0/WrKJK+eMYsna/SxZu58d1Y14BGaVF7K8spXzf17BolllfOG8iUwcZoX3GmP45YtbuOflbZwxfohlb7d/00994SxueHAFNz60im8sPImbzp3Q7UPN0+/t541tR7hz0XSG5kX/QJNMuPYLE5E04D7gQqASWC4ii40xG8KaXQfUGmMmisjVwE+Bj7kxnhvPnUDm0b0055fz9o4j/O3t3fzx9Z2IgEeEKcPz+L/PzafQF5taNmZIDj+/6hQA9tU1sXxnDct31bBiVy13v7ClvV26Rxg9xMf4klzGDPGxvuooy3bW0GZgzBAf150zjoXTh+O3E3u9sOEgL244iEfAm+5hWj8X/872pjF3bBHPrN1PZrqHWy6OzEbc17UOHo/wi6tmcdm9r7v+FBoNZ00saa/ZkAjMKMtnRlk+eW2NcVu5HM5Jw/P4yKmj+NfqfVx68oi4998dV5w6ivcqj/LnN3a2C7rTxw/hs2eP4+Lpwxmal8lTz7/M+tbhPPj2bp5avY8PzRzJDe8bzwOv7mDxmio+OmcUP/zwyR20pKF5mTx8/enc8vh73PXcZnYcbuRHndoANLS08oNnNjCjLL89KCKVcfNRax6wzRizA0BEHgEWAeGCYBFwh/36ceBeERHjQvWKrIw0pg5JY8EC60muORhizd463tlZQ01jgK9eMClmIdCZssJsymaXcfnsMsBalr/9cAM7Djeyw/6/s7qRV7ceZkyxjy+eN5GFM0YwdUReh6eTsyaWcMdl01lfdYwX1h/gpY2HOCPKtAHx4JxJQ3lz+xE+d844RhZG52jrC0PzMnn56wvI7KcoqWRERHjipjN57dVXXbvGDy6fwfXnjO/3p+LbPjCV3Mx0RhRmcfH04SeYGQszPdx28TQ+f+4E/vDqDv721u52X0ZPT/tZGWncc/UsJgzN4VcvbeXZtftPiGprDbXhD4b43SfnRJ1KIxkRtyoGiciVwEJjzOfs7U8B840xXwxrs85uU2lvb7fbVHfq6wbgBoDS0tI5jzzySExjamhoIDc3t/eG/UT4Cku3iMec61raeGZHkCsmeclOT/wfRaJ9zv2BzhmOBQxL9wQZne9h9rDInnHfPdTK+uqu8wJNLErj9BGJYZZ06MvnfN555600xszt6lhizbIbjDEPAA8AzJ071yxYsCCmfioqKoj13GQlXnO+vM899B/6OQ8OuprzZVH2saDXFomFW5+zmzr3PiC8NNIoe1+XbUQkHSjAchoriqIo/YSbgmA5MElExomIF7gaWNypzWLgM/brK4GX3fAPKIqiKN3jmmnIGNMqIl8EnscKH/2zMWa9iNwJrDDGLAb+BDwoItuAGixhoSiKovQjrvoIjDFLgCWd9t0e9roZ+KibY1AURVF6RuPyFEVRBjkqCBRFUQY5KggURVEGOSoIFEVRBjmurSx2CxE5DOzutWHXlADVvbZKLXTOgwOd8+CgL3MeY4wZ2tWBpBMEfUFEVnS3xDpV0TkPDnTOgwO35qymIUVRlEGOCgJFUZRBzmATBA8M9AAGAJ3z4EDnPDhwZc6DykegKIqinMhg0wgURVGUTqggUBRFGeQMGkEgIgtFZLOIbBORWwd6PPFCRP4sIofsam/OvmIReVFEttr/i+z9IiL32O/BeyJy6sCNPDZEpFxElorIBhFZLyJfsfen8pyzRGSZiKyx5/w9e/84EXnHntujdrp3RCTT3t5mHx87oBPoAyKSJiLvisjT9nZKz1lEdonIWhFZLSIr7H2uf7cHhSAQkTTgPuASYBpwjYhMG9hRxY2/Ags77bsV+I8xZhLwH3sbrPlPsv9uAH7bT2OMJ63A140x04DTgS/Yn2Uqz7kFeL8x5hRgFrBQRE4Hfgr80hgzEagFrrPbXwfU2vt/abdLVr4CbAzbHgxzPs8YMytsvYD7321jTMr/AWcAz4dtfwv41kCPK47zGwusC9veDIywX48ANtuvfw9c01W7ZP0D/gVcOFjmDPiAVcB8rBWm6fb+9u84Vg2QM+zX6XY7GeixxzDXUfaN7/3A04AMgjnvAko67XP9uz0oNAKgDNgbtl1p70tVSo0x++3XB4BS+3VKvQ+2+j8beIcUn7NtIlkNHAJeBLYDdcaYVrtJ+Lza52wfPwoM6dcBx4dfAd8A2uztIaT+nA3wgoisFJEb7H2uf7eToni9EjvGGCMiKRcjLCK5wBPAV40xx0Sk/VgqztkYEwJmiUgh8E9gysCOyF1E5IPAIWPMShFZMMDD6U/ONsbsE5FhwIsisin8oFvf7cGiEewDysO2R9n7UpWDIjICwP5/yN6fEu+DiGRgCYH/M8Y8ae9O6Tk7GGPqgKVYZpFCEXEe5sLn1T5n+3gBcKR/R9pnzgIuE5FdwCNY5qFfk9pzxhizz/5/CEvgz6MfvtuDRRAsBybZEQderNrIiwd4TG6yGPiM/fozWHZ0Z/+n7WiD04GjYSpnUiDWo/+fgI3GmF+EHUrlOQ+1NQFEJBvLJ7IRSyBcaTfrPGfnvbgSeNnYRuRkwRjzLWPMKGPMWKzf68vGmE+QwnMWkRwRyXNeAxcB6+iP7/ZAO0f60QlzKbAFy7Z620CPJ47zehjYDwSxbITXYdlG/wNsBV4Ciu22ghU9tR1YC8wd6PHHMN+zseyo7wGr7b9LU3zOM4F37TmvA263948HlgHbgH8Amfb+LHt7m318/EDPoY/zXwA8nepztue2xv5b79yn+uO7rSkmFEVRBjmDxTSkKIqidIMKAkVRlEGOCgJFUZRBjgoCRVGUQY4KAkVRlEGOCgIlaRGRkJ2l0fkbO9Bjigcicq2IHBaRP9rbC5zsm2Ft/ioiV3bdA4jIz0TkgIj8t9vjVZIfTTGhJDNNxphZXR2wF56JMaatq+NJwKPGmC/GerIx5hYRaYzngJTURTUCJWUQkbFi1Zz4G9bCq3IRuUVEltv52r8X1vY2EdkiIq+LyMPOk7OIVIjIXPt1iZ3iwEn69rOwvj5v719gn/O4iGwSkf+zhRAicpqIvClWHYFlIpInIq+KyKywcbwuIqf0Yc5zwzSitamWY0npH1QjUJKZbDsjJ8BO4GtYudk/Y4x5W0QusrfnYa3CXCwi7wMasdIWzML6DawCVvZyreuwlvCfJiKZwBsi8oJ9bDYwHagC3gDOEpFlwKPAx4wxy0UkH2jCSo9xLfBVEZkMZBlj1kQw13PC5gowGmu17Qp7HojIz4DnIuhLUTqggkBJZjqYhmwfwW5jzNv2rovsv3ft7VwswZAH/NMY47fPiyTv1EXAzDC7fIHdVwBYZoyptPtajVUf4iiw3xizHMAYc8w+/g/gOyJyC/BZrMJCkfCaMeaDYXPtcJ6IfAw41R6nokSFCgIl1Qi3iwvwY2PM78MbiMhXezi/leMm06xOfX3JGPN8p74WYFUQcwjRw+/KGOMXkReBRcBVwJwexhIRIjIDuAN4n7HSVStKVKiPQEllngc+K1btAkSkzM7z/ipwuYhk29kePxR2zi6O35yv7NTXTXYKbERksp0hsjs2AyNE5DS7fV5Y+uQ/AvcAy40xtX2ZoJ2V9GHg08aYw33pSxm8qEagpCzGmBdEZCrwlu2/bQA+aYxZJSKPYmV5PISVptzhbuAxsapDPRO2/49YJp9VtjP4MHB5D9cO2Oaa39ipo5uAC4AGYxVbOQb8JQ7TXASMAf5gz5HuIqkUpTs0+6gy6BGRO7Bu0Hf30/VGAhXAlK7CW0XkWqyUwjGHj9r93EE/zktJXtQ0pCj9iIh8GqvG8m09rHFoAi5xFpTFeJ2fAZ+ko89EUbpENQJFUZRBjmoEiqIogxwVBIqiKIMcFQSKoiiDHBUEiqIogxwVBIqiKIOc/w+g6XqzUOhDVQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "f, Cxy = signal.coherence(x, y, fs, window='hamming',nperseg=100, noverlap=80)\n",
    "# noverlap段之间重叠的点数 若None,noverlap = nperseg // 2.默认为None.\n",
    "\n",
    "plt.plot(f, Cxy)\n",
    "plt.title('Magnitude-Squared Coherence')\n",
    "plt.xlabel('Frequency [Hz]')\n",
    "plt.ylabel('Coherence')\n",
    "plt.grid()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5d2fa8ad",
   "metadata": {},
   "source": [
    "当相干性很小时，忽略交叉频谱。绘制交叉频谱的相位，并指出在两个时间点间具有显著相干性的频率。标记已知的正弦分量间相位滞后。在 100 Hz 和 200 Hz 时，根据交叉频谱估计的相位滞后接近真实值。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "aa2d3602",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEWCAYAAABrDZDcAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAtF0lEQVR4nO3deZxkZX3v8c+3l+keYBZhYEQYGRASM0aDMgIqVxqCCbiASdCAxhVFr6ImiNu9XoLEGMXXdbmRJCIajRuC64RggAiNUQMOIMquAxkYBmQbmJ6tl+r63T/OOT01Pb1UV53TXafr+3695jVVp06dep7q7vrV83s2RQRmZta+Oua6AGZmNrccCMzM2pwDgZlZm3MgMDNrcw4EZmZtzoHAzKzNORCYtQlJ/ZLeMtflsNbjQGC5kPQaSTdK2irpIUk/lHTMHJbnGEk/k7RZ0iZJP5X0/IJfc72kE4p8jTrLsCP9OTws6cuS9prLMlnrcyCwpkk6G/gM8DFgOfB04B+AUyY5v6vg8iwGLgf+HtgbOAD4CDBU5OvWUa5C613jFRGxF/A8YDXw4Vl6XSspBwJriqQlwPnAOyPiuxGxLSJGIuJfI+J96TnnSfq2pK9JGgDeKOlpktak39bXSXprzTWPTFsXA+m32k+lx3vTazwu6UlJayUtn6BYvwMQEd+MiNGI2BERV0XEr9LrvDFtIXwubTHcJekPa+sk6Ytpy2ajpI9K6qx5/K2S7pS0RdIdkp4n6askAfBf02/j75e0UlJIOkPS/cA1kvokPTDuPRxrSaTv1WVpPbdIulXS70j6kKRHJG2Q9Ef1/GwiYiPwQ+D3aw4flNZ9i6SrJC2rKcdlkn6bvic/lvSsmsdemtZ1S/qenFPz2Msl3ZL+TH4m6Tn1lM9ahwOBNesFQC/wvWnOOwX4NrAU+DpwCfAA8DTgVOBjko5Pz/0s8NmIWAw8A7g0Pf4GYAmwAtgHeDuwY4LX+jUwKukrkk6S9JQJzjkKuAdYBvw18F1Je6ePfRmoAIcCzwX+CHgLgKRXAecBrwcWAycDj0fE64D7Sb+NR8QFNa91LPB7wB9P8x5lXgF8FXgK8AvgSpK/1QNIgu7n67mIpBXAS9NrZF4DvAnYD1gAnFPz2A+Bw9LHbib5OWW+CLwtIhaRBJZr0td4LvAl4G0kP5PPA2sk9dRZV2sBDgTWrH2AxyKiMs15/xUR34+IKsmH74uAD0TEYETcAlxM8uEKMAIcKmlZRGyNiOtrju8DHJp+078pIgbGv1B67BgggC8Aj6atj9rWwyPAZ9LWy7eAu4GXpee8FPjLtHXzCPBp4LT0eW8BLoiItZFYFxH3TVP389JrTRS0JvKfEXFl+p5eBuwLfDwiRkgC6EpJS6d4/vclPQn8BLiOJGWX+eeI+HValkuBw7MHIuJLEbElIoZIgt0fpC0+SN77VZIWR8QTEXFzevxM4PMRcUP6M/kKSQru6Drrai3AgcCa9TiwrI7894aa208DNkXElppj95F84wU4gyS9c1ea/nl5evyrJN+OL5H0oKQLJHVP9GIRcWdEvDEiDiT5Bvs0kn6MzMbYdcXF+9JzDgK6gYfSVMeTJN9y90vPW0HSkpiJDdOfsouHa27vIAm0ozX3AabqAH5lRCyNiIMi4h3jAtBva25vz64jqVPSxyXdk6bv1qfnZKmjPyMJkPdJuk7SC9LjBwHvzd6r9P1aQfJeWkk4EFiz/ovkG+Arpzmv9kP3QWBvSYtqjj0d2AgQEb+JiNNJPnw/AXxb0p7pt/ePRMQq4IXAy9nZipj8hSPuIkn31ObKD5Ckca//IMmH9hCwLP0wXRoRiyMiy5dvIElXTVfHyY5vA/bI7qR9D/tOV4dZ8BqS9N0JJOm3lelxAaQtoFNIfibfZ2e6bgPwtzXv1dKI2CMivjmbhbfmOBBYUyJiM3AucKGkV0raQ1J3mpu/YJLnbAB+Bvxd2gH8HJJWwNcAJP2FpH3TNNKT6dOqko6T9Oz0w3OAJF1RHX99Sc+U9F5JB6b3VwCnA9fXnLYf8O60rK8iyeFfEREPAVcB/1fSYkkdkp4h6dj0eRcD50g6QolDJR2UPvYwcMg0b9mvgV5JL0tbMx8GWiGfvogkAD5OEqjG0kmSFkh6raQlaXpqgJ3v+xeAt0s6Kn0/9kzrtmj8C1jrciCwpkXE/wXOJvlQe5TkW+JZJN8cJ3M6ybfOB0k6mv86Iv4jfexE4HZJW0k6jk9L0xtPJelwHgDuJMl/f3WCa28h6Qy+QdI2kgBwG/DemnNuIOkYfQz4W+DUiHg8fez1JB2pdwBPpK+5f1rXy9Lzv5G+zvdJhqgC/B3w4TRFUtsJW/tebQbeQRJQNpK0EB6Y6NxZ9i8k6bGNJPW+ftzjrwPWp2mjtwOvBYiIG4G3Ap8jea/WAW+cnSJbXuSNaazdSHoj8JaImLMJb2atxC0CM7M250BgZtbmnBoyM2tzbhGYmbW5ohf/OpFk1EcncHFEfHzc428EPkk6fhz4XERcPNU1ly1bFitXrmyoPNu2bWPPPfds6Lll5Tq3B9e5PTRT55tuuumxiJhwzkphgSAd630h8BKS4XFrJa2JiDvGnfqtiDir3uuuXLmSG2+8saEy9ff309fX19Bzy8p1bg+uc3tops6SJl0KpcjU0JHAuoi4NyKGSdZImXBZYjMzmzuFdRZLOhU4MSKyVRtfBxxV++0/TQ39HckkpF8Df5XOOh1/rTNJFrdi+fLlR1xyySUNlWnr1q3stVd77dHhOrcH17k9NFPn44477qaIWD3RY7O1UcZk/hX4ZkQMSXob8BXg+PEnRcRFwEUAq1evjkabRm5KtgfXuT24zvkpMjW0kWQVwsyB7OwUBiAiHk+XvIVkyv0RBZbHzMwmUGQgWAscJulgSQtI1nNfU3uCpP1r7p5Msn6MmZnNosJSQxFRkXQWyfrxncCXIuJ2SecDN0bEGpLVH08m2Q1qE16sysxs1hXaRxARVwBXjDt2bs3tDwEfKrIMZmY2tbnuLDbLxfbhCnc+tIU7HtzMwGOj9M11gcxKxIHASum/H9vGj+58mNs2bua2Bwe499GtVGtGQvc+9b8545iD566AZiXiQGCl9N5Lb+Hm+59k/yW9POtpS3jZs/fn9w9Ywu8uX8TZ//Jj/ubyO9i0bYhz/uh32XVHSjMbz4HASmnLYIUTn/VU/ul1u484fsfhPVz9xDIuvPYeNm0b5qOvfDadHQ4GZpNxILBSGqpUWbigc8LHOiQ+9ie/zz57LuBz167jiW0jfOa0w+ntnvh8s3bnZaitlIYqo/R0Tf7rK4lz/vh3Offlq/j323/Lm/55LVsGR2axhGbl4UBgpTRUqU4ZCDJvPuZgPvPnh7N2/Sbe9c1fzELJzMrHqSErpaGRKj11pnpe+dwD+Mm6x/jP3zxacKnMysktAiudiJg2NTTengs6GapUCyyVWXk5EFjpVKpBNZhRIOjp7mRoxIHAbCIOBFY62Tf7nq76RwH1dnUwWBmlqP03zMrMgcBKZ2hkFICe7pm1CCJgZNSBwGw8BwIrnaxFsKBzBoEgTSMNVUYLKZNZmTkQWOmMpYZm0iIYCwTuJzAbz4HASif7Vj+TPoLsXAcCs905EFjpDI91Fs+kjyBtEYw4NWQ2ngOBlU4jo4acGjKbnAOBlU42H2BmfQRODZlNxoHASmdnH8HMO4sHnRoy240DgZVOQ6mhbrcIzCbjQGCl00yLwJ3FZrtzILDSaaSPoLfbncVmk3EgsNJpbNSQU0Nmk3EgsNJpKjXkJSbMduNAYKUzlhqaUSDo3OW5ZraTA4GVzlClSmeH6JrJonPuIzCblAOBlc5MdyeDnSuVOjVktjsHAiudoUqVBTMMBB0dYkFnB4NODZntxoHASmdopDrjFgEk6SG3CMx250BgpTM8Wp3R0NFMT5c3sDebiAOBlU4jfQSQjDLyqCGz3TkQWOkMjVRnNKs449SQ2cQKDQSSTpR0t6R1kj44xXl/JikkrS6yPDY/DFWcGjLLU2GBQFIncCFwErAKOF3SqgnOWwS8B7ihqLLY/NJUasiBwGw3RbYIjgTWRcS9ETEMXAKcMsF5fwN8AhgssCw2jyQtgkb7CJwaMhuvq8BrHwBsqLn/AHBU7QmSngesiIh/k/S+yS4k6UzgTIDly5fT39/fUIG2bt3a8HPLaj7WedOT2+ke6Zi0XpPVeduWQbaPxLx7P2B+/pyn4zrnp8hAMCVJHcCngDdOd25EXARcBLB69ero6+tr6DX7+/tp9LllNR/r3LX2Wg7cfyl9fc+d8PHJ6vz1+29kw6bt9PW9uOASzr75+HOejuucnyJTQxuBFTX3D0yPZRYBvw/0S1oPHA2scYexTafR1FBvdyfD7iMw202RgWAtcJikgyUtAE4D1mQPRsTmiFgWESsjYiVwPXByRNxYYJlsHmh81JA7i80mUlggiIgKcBZwJXAncGlE3C7pfEknF/W6Nv8NjTQzasidxWbjFdpHEBFXAFeMO3buJOf2FVkWmz+GKg1OKOvq9Mxiswl4ZrGVSmW0SqUajaWGup0aMpuIA4GVyvBo8kE+02WoIUkNDY9WqVYj72KZlZoDgZVKI9tUZrJWRBZMzCzhQGClkn2INzpqCGDQs4vNduFAYKXSTIugtzvdwN79BGa7cCCwUsmGfzY2aijdt9gjh8x24UBgpZJ9m2901FByDaeGzGo5EFipjLUImugsdmrIbFcOBFYqzY0acovAbCIOBFYqY6mh7sZHDbmPwGxXDgRWKk2lhjxqyGxCDgRWKjs7i50aMsuLA4GVylgfQROpoUGnhsx24UBgpdJMamjnhDK3CMxqORBYqeSTGnKLwKyWA4GVSnMTytIWgVNDZruoOxBI2lPSzP/6zHI0lC4Y192pGT/XncVmE5s0EEjqkPQaSf8m6RHgLuAhSXdI+qSkQ2evmGaJbON6aeaBoKtDdMipIbPxpmoRXAs8A/gQ8NSIWBER+wHHkGw0/wlJfzELZTQbkwWCRkhKtqt0IDDbxVR7Fp8QESPjD0bEJuA7wHckdRdWMrMJJPsVN56h7OnuGEsvmVliqkCwaKrmd0RsmihQmBVpqDLacIsAkn4CtwjMdjVVILgJCEDA04En0ttLgfuBg4sunNl4zaSGIBlt5B3KzHY16V9URBwcEYcA/wG8IiKWRcQ+wMuBq2argGa1hkaqDQ0dzfR2u0VgNl49X62OjogrsjsR8UPghcUVyWxyQ5XRhnYny7iz2Gx3U6WGMg9K+jDwtfT+a4EHiyuS2eSaTw11eB6B2Tj1/EWdDuwLfC/9t196zGzWJYGg2VFDbhGY1Zq2RZAOF33PLJTFbFpDI6P0LOpp+Pk9XZ0M7KjkWCKz8ps2EEjaF3g/8CygNzseEccXWC6zCQ03O4/AqSGz3dSTGvo6yfISBwMfAdYDawssk9mk8ukjcGrIrFY9f1H7RMQXgZGIuC4i3gy4NWBzovkJZZ5HYDZePaOGstnDD0l6GcmIob2LK5LZ5IZGqixoJhB4HoHZbuoJBB+VtAR4L/D3wGLgrwotldkkmh011Nvd6VFDZuNM+dUq3X/gsIjYHBG3RcRxEXFERKyp5+KSTpR0t6R1kj44weNvl3SrpFsk/UTSqgbrYW2gWg2GR/OZRxAROZbMrNym/IuKiFEanDOQBpELgZOAVcDpE3zQfyMinh0RhwMXAJ9q5LWsPQyPZhvXNxcIqgGVqgOBWaae1NBPJX0O+BawLTsYETdP87wjgXURcS+ApEuAU4A7aq4xUHP+niSL3JlNqJltKjPZc4cqVbo7vVOrGdQXCA5P/z+/5lgw/cihA4ANNfcfAI4af5KkdwJnAwsmu6akM4EzAZYvX05/f38dxd7d1q1bG35uWc2nOj85lASC++5dR3/lvknPm6rO992XjH245rr/ZPGCme9y1qrm08+5Xq5zjiKikH/AqcDFNfdfB3xuivNfA3xluuseccQR0ahrr7224eeW1Xyq8/2Pb4uDPnB5XLr2/inPm6rOl/z8vjjoA5fHxie251y6uTWffs71cp1nBrgxJvlcLbJtvBFYUXP/wPTYZC4BXllgeazkxlJDTc0s3pkaMrNEkYFgLXCYpIMlLQBOA3YZbSTpsJq7LwN+U2B5rOSypSGaHTUEeFKZWY16+ggaEhEVSWcBVwKdwJci4nZJ55M0UdYAZ0k6gWTS2hPAG4oqj5Xfzs7i5iaU1V7LzKYIBJKOj4hrJP3pRI9HxHenu3gkG9pcMe7YuTW3vaqp1S2bCNbUhLIsNeQWgdmYqVoExwLXAK+Y4LEApg0EZnkaSw01M4/ALQKz3UwaCCLir9P/3zR7xTGbXC6pIXcWm+1m2r8oSftI+n+SbpZ0k6TPStpnNgpnViufCWVZi8CpIbNMPV+tLgEeBf6MZG7AoySzjM1mVZbXz6VF4IXnzMbUM2po/4j4m5r7H5X050UVyGwyHjVkVox6/qKuknSapI7036tJhoSazSqnhsyKUU8geCvwDWAIGCZJFb1N0hZJA1M+0yxHuYwaSoPIoFNDZmOmTQ1FxKKJjktaEBHD+RfJbGLDaYtgQROrhrpFYLa7Gc0slnQQyf4CJwH7kyw1bTYrhipVFnR20NHR+KqhHR1iQae3qzSrNWUgkNRDMrHspcAxJMtKXwm8JyLWF146sxpDI83tTpbp6erwqCGzGlMtMXEFsA9wHclice+LiJHJzjcr2lBltKn+gUyygb1TQ2aZqVoEb42IXZaNlvQUYEVE/KrYYpntrtmN6zM9XZ1ODZnVmPTrVRYEJPVLWixpb+Bm4AuSvLewzbokEOSUGnIgMBtTz1/Vkkj2Fv5T4F8i4ijghGKLZba7oZFRFuQQCBZ0dXj1UbMa9fxVdUnaH3g1cHnB5TGb1FCl2tTuZJme7k4G3SIwG1NPIDifZKTQPRGxVtIheCcxmwNDldEcRw25RWCWqWdC2WXAZTX37yVZgM5sVg1VquzV0/ymer3dnQzs8AA4s0w9y1D/jqQfSbotvf8cSR8uvmjWjNFq8N2bH+AH64Z5ZGBwrouTi2QeQR6jhtxZbFarnnb2F4APkewrTDp09LQiC2WNiwiuvuNhTvrsjzn70l/yvXUj/I8LruW8Nbfz0OYdc128puQ2j6DL8wjMatXTzt4jIn4u7TKtv1JQeawJa9dv4hM/vIsb73uCQ5btyT+89nkM3H8nN+/Yh69dfx/fuOF+Xv38A/mffYdywNKFc13cGRuqVOlpYp2hTE9Xp2cWm9WoJxA8JukZJPsUI+lU4KFCS2UzsmHTds5bczs/uusR9lvUw8f+5Nm8avWBdHd20P/43Vzwsj/gXccfxj9edw/fWruBb63dwJkvPoT3/fEz57roM5KMGsprZrEDgVmmnkDwTuAi4JmSNgL/DfxFoaWyGfnklXfz03se431//Lu8+UUHs3DB7nn0FXvvwcf+5Nmcddyh/K/v3cqF197Du44/jN4chmPOlqGR0Rz7CJwaMsvUM2roXuAESXsCHRGxpfhi2Uw8smWQ5xywlHced+i05z5t6UL+8Jn70X/3o2wZrJQqEAyP5jWz2EtMmNWaNhBIWgq8HlhJMrkMgIh4d5EFs/oN7KjwtKW9dZ+/eGF38rzBEfZd1FNUsXIVEbkuMTFcqVKtRlNLWpvNF/Wkhq4ArgduBfw1qgUNDI7wzIUT7h80oSwQbC7RWPqR0SCCnGYWJ8FkeLRKb0d5WkRmRaknEPRGxNmFl8QatnnHCIt7u+s+Pzu3TJOqxrapzKFF0Jv2MwyNVEuVGjMrSj1/VV+V9FZJ+0vaO/tXeMmsLtVqsHWoMvYtvx5LFibxf2CwPKOAd25cn8+ooeSa7jA2g/paBMPAJ4H/TTqENP3/kKIKZfXbMlQhApbMIBCUs0WQBYJ89iOovaZZu6snELwXODQiHiu6MDZz2Yf54t761+ApYx9BtkhcXjOLwS0Cs0w9f1XrgO1FF8Qak32YzyQ11NvdyYKuDgYGSxQI8kwNpdcY9OxiM6C+FsE24BZJ1wJD2UEPH20N2Yf5TFJDkKSHBnaUsY8gn/0Iaq9p1u7qCQTfT/9ZC8o+zGcyagiSDuNStQhG8hs15NSQ2a7qmVn8lUYvLulE4LNAJ3BxRHx83ONnA28hWcTuUeDNEXFfo6/Xjsb6CBbObJ3+xQu7y9lZnGcfgVNDZsAUfQSSLk3/v1XSr8b/m+7CkjqBC4GTgFXA6ZJWjTvtF8DqiHgO8G3ggkYr0q6yb/Uz6SOALDVUwkCQ66ghtwjMYOoWwXvS/1/e4LWPBNalaxUh6RLgFOCO7ISIuLbm/OvxYnYzNrBjhA7BXgtm3iK4f1N5xgBkH9p5bF7fOzaPwC0CM5giEEREttT0OyLiA7WPSfoE8IHdn7WLA4ANNfcfAI6a4vwzgB9O9ICkM4EzAZYvX05/f/80Lz2xrVu3NvzcVnX7uiEWdsGPf3zdhI9PVudtm4Z4bKBSmvfjlxuT1sstN63lwT2mDgbT/Zwf35EEgF/edgdLnpwf22/Px9/t6bjO+anna+RL2P1D/6QJjjVM0l8Aq4FjJ3o8Ii4iWQqb1atXR19fX0Ov09/fT6PPbVU/ePgW9t66adJ6TVbnnw/exY833suxxx7LuE2HWtKDN9wPt97Ki1/0Qp66ZOoF9qb7OT+2dQiu+w8OfsZh9L1gZb4FnSPz8Xd7Oq5zfiYNBJL+J/AO4JBxfQKLgJ/Wce2NwIqa+wemx8a/zgkks5aPjYih8Y/b1AZ2jMx46CgkqaFKNdgxMsoeM0wrzYXhHNca2jlqyKkhM5i6RfANklTN3wEfrDm+JSI21XHttcBhkg4mCQCnAa+pPUHSc4HPAydGxCMzKbglZrrgXGbnMhOVUgSCfEcNeR6BWa2pPgFGI2I9cPpkJ0jaKyK2TvRYRFQknQVcSTJ89EsRcbuk84EbI2INyRpGewGXpemJ+yPi5Maq0p4GBkc4ZNleM37ekpplJqZLtbSC7EN7QQ57Fnd3Cmnn3ASzdjdVIPiBpFuAHwA3RcQ2AEmHAMcBrwa+QDLsc0IRcQXJfga1x86tuX1CwyU3IPlG31hqKFuBtBxDSIcqo3R1iK4cAoEkero6GHSLwAyYetTQH0p6KfA24EWSnkIy8etu4N+AN0TEb2enmDaZzTtGZjyZDMq3AunQSD67k2V6ujrdIjBLTfkJMtE3emsdw5UqO0ZGG+ojWFKzXWUZDFWquexOlkk2sHeLwAzqW33UWtSWbMG5PRobNQSweXtZAsFori2C3m5vYG+WcSAosbElqBtoESzqLdcuZXltXJ9JWgRODZmBA0GpZR/ijfQRdHd2sMeCzpL1EeSYGuru8KJzZqlpP0Em2Z94S0SU4xNkHss+xBsZNZQ9rzx9BKO5zCHI9HQ5NWSWqecv62aSJaJ/Dfwmvb1e0s2SjiiycDa1sZVHG0gNZc8ry3aVTg2ZFaeev6yrgZdGxLKI2IdknaHLSZaf+IciC2dTa2SbylqLF3aVZpeyJBB41JBZEeoJBEdHxJXZnYi4CnhBRFwP9BRWMptW9iHeaGpocW+5UkN5LEGd6enqZNDzCMyA+gLBQ5I+IOmg9N/7gYfTjWf8lWoODQyOsKCzo+GUSan6CPKeUNbtFoFZpp6/rNeQrBz6/fTf09NjnSTLTNgcyWYVN7qM9OKF3SWaR1BAH4FHDZkB9e1Z/BjwrkkeXpdvcWwmBhpceTSzuLeLLUMVqtWgo6O19yQYzrmPIJlQ5tSQGdQ3fHRf4P3As4CxZSoj4vgCy2V1GBisNNxRDEmLIAK2DleaCiizIf/ho04NmWXq+cv6OnAXcDDwEWA9yV4DNseS1FBzgQDKscxE/qkhzyMwy9Tzl7VPRHwRGImI6yLizYBbAy1gy44RFvc2vqnM2AqkJegwLmL46Gg1qIw6GJjV8ymSfUo8JOllwIPARLONbZYNDDa2TWVmbE+CFp9LUBmtMlqN3EcNQRJg8tjjwKzM6gkEH5W0BHgv8PfAYuAviyyUTS8imk4NlWUp6jy3qczUble5p2fDWJurZ9TQ5enNzSQ7kyHpLwssk9VhcKTKyGg0OWpo53aVrWwsEOScGgI8qcyMxlcfPTvXUtiMZd/im0sNlWOXsmyYZ1GpIbN21+hfVmsPOm8DO9cZaryzeFFPF1Lr70mQTfzKMzXUO5YacovArNG/rMi1FDZjA01sSpPp6BCLerpK0CIoIDWUtQg8u9hs8j4CSVuY+ANfwMLCSmR1ySM1BEl6qPUDQQGpoZrOYrN2N2kgiIhFs1kQm5lsyGczo4agHCuQFtlZ7NSQmbeqLK2d+xU33kcA6QqkLT6PIEvf5L0Mde21zdqZA0FJDTS5KU1m8cKuErQIPGrIrEgOBCU1MDjCHgs66W5yVmwZtqssZkKZU0NmGQeCktrc5BLUmTJ0Fg8X0keQXGvQqSEzB4KyGthRaXrEECR9BNuGR1t68bViRg25RWCWcSAoqYHBkaYmk2WyzuZWnlS2c9RQjhPKuj181CzjQFBSeaaGoLWXmdg5szi/1FA2AsmjhswcCEqr2SWoM2VYgbSI1FBnh+julFNDZjgQlNbAjua2qczsbBG0dmqoQ9CV877K3qXMLFFoIJB0oqS7Ja2T9MEJHn+xpJslVSSdWmRZ5pNqNZI+giYnk0E5lqLOdieT8g4EHW4RmFFgIJDUCVwInASsAk6XtGrcafcDbwS+UVQ55qOtwxUimp9MBjW7lLVyamgk343rMz1dHe4jMKO+HcoadSSwLiLuBZB0CXAKcEd2QkSsTx/zX+MM5DWrGGr6CFq+RVBAIOjuZNCpIbNCA8EBwIaa+w8ARzVyIUlnAmcCLF++nP7+/oYKtHXr1oaf20ruG0jSGRvuuZv+rfdMee50dY4IOgW33n0P/bFh0vPm0n0bB6mOVOv+2dX7cx4Z3M7Ghwbnxe/EfPndngnXOT9FBoLcRMRFwEUAq1evjr6+voau09/fT6PPbSX/dc/j8LPrecHqw3nhM5ZNeW49dV7yk6tZut9T6et7do6lzM+lG29iaWUrfX3H1nV+vT/nfW77CYv2WEBf35FNlnDuzZff7ZlwnfNTZGfxRmBFzf0D02PWpLz2Isgs7u1q7VFDI9Vi+gi6O91ZbEaxgWAtcJikgyUtAE4D1hT4em0jj93Jai1Z2Np7EgxVqixocnG9iSSjhtxHYFZYIIiICnAWcCVwJ3BpRNwu6XxJJwNIer6kB4BXAZ+XdHtR5ZlPNufYWZxdp7U7i0dzXXAu09PV6VFDZhTcRxARVwBXjDt2bs3ttSQpI5uBgcEKUrL5fB4W93bz4JM7crlWEYYqVfbcM/9f1Z5uzyMwA88sLqWBHSMs6umiI6eZtosXdrf0onPDRQ0fdWrIDHAgKKWBHSO5pYUg3aWspVND1eJSQw4EZg4EZZTXgnOZxb3dDFWqDI60ZppkaGS0sBZBq9bZbDY5EJTQwI5KbiOGoGbhuRYdOTRUKWr4qFNDZuBAUEqbd+SzKU1mSYuvQFpkami4UiUicr+2WZk4EJRQ/qmh1l54Lhk+mv+vam93tl2lWwXW3hwISmggp93JMllqqBWXoh6tBiOjUViLABwIzBwISmZktMq24dF8Rw31tu4KpMPZfsUFLUMN3sDezIGgZLak4/3zTA3t3K6y9foIitimMtPjfYvNAAeC0tm5F0F+ncWLsj6CFmwRZGmbQlJD3U4NmYEDQelsznnBOYDe7k56ujpaMxCMZIHAqSGzojgQlEw2sifPPgJo3RVIx1JDBfYRDDo1ZG3OgaBksrH+efYRQLYCaSv2ESQf0sUsQ52lhtwisPbmQFAyYy2CHFNDyfW6WrxFUEQfgecRmIEDQelsLqCzOLled0vOIyiyj6A3axE4NWRtzoGgZAZ2jNDdKRbm/A15SYtuTjM0WmBncbc7i83AgaB0BgaTWcVSPnsRZBb3tuaeBDtbBEXMLHZqyAwcCEpn845K7iOGIEk1bd4x0nILsBU7asjzCMzAgaB08t6UJrO4t5vRarB9uLXSJDsnlBWYGvKeBNbmHAhKJkkN5b9/75IW3ZOg0JnFTg2ZAQ4EpbO5qBZBi+5JkH1bLyI1lM1NcIvA2p0DQckM7KjkPpkMds5LaLUhpEWmhiR5A3szHAhKJxs1lLedu5S1ZiAoYmYxJOssORBYu3MgKJHBkVGGK9XcJ5PBzglqrddHkOxOlvdw2UzSInBqyNqbA0GJZN/Wi0wNtVyLYKRaSFoo09Pd4ZnF1vYcCEqkqHWGYOeeBJtbrbO4Ui1knaFMT5dTQ2YOBCWyc52h/ANBV2cHe/W03sJzRW1cn3FqyMyBoFSKWoI6s7i3q/VSQ5UqCwoPBG4RWHtzICiRnamh/DuLId2ToNVaBCPVQiaTZXq6Ohn0PAJrcw4EJVJkagiSvofWm0dQcGqo2y0CMweCEhkoYL/iWq24S9lwpeBRQ10eNWTmQFAiA4MVFnZ3FpYzX7ywFTuLix01lEwoc2rI2luhgUDSiZLulrRO0gcneLxH0rfSx2+QtLLI8pRdsvJoMf0D0KqpoVloETg1ZG2usL8wSZ3AhcBJwCrgdEmrxp12BvBERBwKfBr4RFHlmQ827yhmeYnMkoXdbB2qUK22zp4ExQ8f9TwCs+K+XsKRwLqIuBdA0iXAKcAdNeecApyX3v428DlJigJ2R7l07QY+85Pt7HnzdXlfetZsfHIHv7f/4sKuv3hhNxHwkk9fR0dBSzrM1P2Pb+fwFUsLu35PVwdPbB/mJZ8q7+8FwLbt5f7dbkQ71vmE/Sv0FXDdIgPBAcCGmvsPAEdNdk5EVCRtBvYBHqs9SdKZwJkAy5cvp7+/f8aFeeDhCst7qnRqx4yf2yqWPAWOXLp9RvXfunVr3ecv2l7l6P07qVQHGytgAZ67XweHdT5eWJ2fNjrK85d3Uo3y/l4A7FXy3+1GtGOdOyqVhj7/plNkIMhNRFwEXASwevXq6Ovrm/E1+oDn9ffTyHPLrH+GdX51cUWZNTOpcx/w5iILM0tm+nOeD1zn/BTZWbwRWFFz/8D02ITnSOoClgCPF1gmMzMbp8hAsBY4TNLBkhYApwFrxp2zBnhDevtU4Joi+gfMzGxyhaWG0pz/WcCVQCfwpYi4XdL5wI0RsQb4IvBVSeuATSTBwszMZlGhfQQRcQVwxbhj59bcHgReVWQZzMxsap5ZbGbW5hwIzMzanAOBmVmbcyAwM2tzKttoTUmPAvc1+PRljJu13AZc5/bgOreHZup8UETsO9EDpQsEzZB0Y0SsnutyzCbXuT24zu2hqDo7NWRm1uYcCMzM2ly7BYKL5roAc8B1bg+uc3sopM5t1UdgZma7a7cWgZmZjeNAYGbW5tomEEg6UdLdktZJ+uBclycvkr4k6RFJt9Uc21vS1ZJ+k/7/lPS4JP2/9D34laTnzV3JGyNphaRrJd0h6XZJ70mPz+c690r6uaRfpnX+SHr8YEk3pHX7VrrcO5J60vvr0sdXzmkFmiCpU9IvJF2e3p/XdZa0XtKtkm6RdGN6rPDf7bYIBJI6gQuBk4BVwOmSVs1tqXLzZeDEccc+CPwoIg4DfpTeh6T+h6X/zgT+cZbKmKcK8N6IWAUcDbwz/VnO5zoPAcdHxB8AhwMnSjoa+ATw6Yg4FHgCOCM9/wzgifT4p9Pzyuo9wJ0199uhzsdFxOE18wWK/92OiHn/D3gBcGXN/Q8BH5rrcuVYv5XAbTX37wb2T2/vD9yd3v48cPpE55X1H/AD4CXtUmdgD+Bmkv2/HwO60uNjv+Mke4C8IL3dlZ6nuS57A3U9MP3gOx64HFAb1Hk9sGzcscJ/t9uiRQAcAGyouf9Aemy+Wh4RD6W3fwssT2/Pq/chbf4/F7iBeV7nNEVyC/AIcDVwD/BkRFTSU2rrNVbn9PHNwD6zWuB8fAZ4P1BN7+/D/K9zAFdJuknSmemxwn+3S7F5vTUuIkLSvBsjLGkv4DvAX0bEgKSxx+ZjnSNiFDhc0lLge8Az57ZExZL0cuCRiLhJUt8cF2c2HRMRGyXtB1wt6a7aB4v63W6XFsFGYEXN/QPTY/PVw5L2B0j/fyQ9Pi/eB0ndJEHg6xHx3fTwvK5zJiKeBK4lSYsslZR9maut11id08eXAI/Pbkmb9iLgZEnrgUtI0kOfZX7XmYjYmP7/CEnAP5JZ+N1ul0CwFjgsHXGwgGRv5DVzXKYirQHekN5+A0kePTv++nS0wdHA5pomZyko+er/ReDOiPhUzUPzuc77pi0BJC0k6RO5kyQgnJqeNr7O2XtxKnBNpEnksoiID0XEgRGxkuTv9ZqIeC3zuM6S9pS0KLsN/BFwG7Pxuz3XnSOz2AnzUuDXJLnV/z3X5cmxXt8EHgJGSHKEZ5DkRn8E/Ab4D2Dv9FyRjJ66B7gVWD3X5W+gvseQ5FF/BdyS/nvpPK/zc4BfpHW+DTg3PX4I8HNgHXAZ0JMe703vr0sfP2Su69Bk/fuAy+d7ndO6/TL9d3v2OTUbv9teYsLMrM21S2rIzMwm4UBgZtbmHAjMzNqcA4GZWZtzIDAza3MOBFZqkt4t6U5JX5/rsuRFUr+SlXJPTu9/WdKp487ZOsXzF6arVw5LWlZ0ea38vMSEld07gBMi4oHag5K6YueaNGX02oi4sZEnRsQOkuUo1udbJJuv3CKw0pL0TySTcH4o6a8knSfpq5J+Cnw1nZH7HUlr038vSp+3j6SrlKztf7Gk+yQtk7RSu+7rcI6k89Lbz5D07+liYP8p6Znp8S+na8L/TNK9td/cJX0gXVv+l5I+nl7j5prHD6u93+B7cH767f8WSRsl/XMz17P25EBgpRURbwceJFm//dPp4VUkLYTTSdam+XREPB/4M+Di9Jy/Bn4SEc8iWc/l6XW83EXAuyLiCOAc4B9qHtufZMbzy4GPA0g6CTgFOCqSfQQuiIh7gM2SDk+f9yag3g/uT9Z84N9S8x6cGxGHk8y+3QR8rs7rmY1xasjmmzVpagTgBGBVzcqki9NVS18M/ClARPybpCemumD6nBcCl9Vcq6fmlO9HRBW4Q1K2RPAJwD9HxPb0dTalxy8G3iTpbODPSRYVq8f7IuLbNWXaWnNbwNeAT0XETXVez2yMA4HNN9tqbncAR0fEYO0JtUtWj1Nh11Zyb811nky/eU9kqPby05TvOyQtkmuAmyIijxUyzwMeiAinhawhTg3ZfHYV8K7sTk1K5sfAa9JjJwFPSY8/DOyX9iH0kKR6iIgB4L8lvSp9jiT9wTSvfTXJN/890ufsnV5rkGQ3rX+k/rTQpCS9gqT18e5mr2Xty4HA5rN3A6uVbOx9B/D29PhHgBdLup0kRXQ/QESMAOeTrF55NVC7KchrgTMkZStDnjLVC0fEv5MsE3xjmtM/p+bhr5PsunVVU7VLnE2yK9XP0/6D83O4prUZrz5qbS8dZrk6Ih6bpdc7B1gSEf9nksf7gXMaHT5ac531zGK9rLzcIjCbRZK+B7yeZETTZDYBX84mlDXwGgvTVkg3O/f7NZuUWwRmZm3OLQIzszbnQGBm1uYcCMzM2pwDgZlZm3MgMDNrc/8fJEcODIOMqmgAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "f, Pxy = signal.csd(x, y, fs,window='hamming', nperseg=100, noverlap=80)\n",
    "for i in range (len(Pxy)):\n",
    "    if Cxy[i]<0.27:\n",
    "        Pxy[i]=0\n",
    "    else:\n",
    "        Pxy[i]=Pxy[i]\n",
    "\n",
    "plt.plot(f, abs(np.angle(Pxy)/np.pi))\n",
    "plt.title('Cross Spectrum Phase')\n",
    "plt.xlabel('frequency [Hz]')\n",
    "plt.ylabel('Lag (times\\pi rad)')\n",
    "plt.grid()\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
